Metadata

title: "E-TRAINEE Tutorial - Sentinel-2 data STAC access and time series processing"
description: "This is a tutorial within the second theme of Module 1 of the E-TRAINEE course."
lastUpdate: 2023-07-13
authors: Andreas Mayr

Sentinel-2 data STAC access and time series processing¶

This Notebook illustrates how to query and download Sentinel-2 satellite imagery from Amazon Web Service (AWS) using a SpatioTemporal Asset Catalog (STAC). Furthermore, it shows how to work with the data using xarray (which was introduced in a tutorial in theme 1), and how to perform tasks like subsetting, cloud masking, spectral index calculation, temporal aggregation, and trend analysis. A few simple processing steps let you explore the development of the Normalized Difference Vegetation Index (NDVI) in spring and summer, including e.g. the seasonal greening of mountain grasslands.

The Notebook heavily relies on the Python packages pystac-client, stackstac, and xarray. These packages are contained in the requirements file provided for the course. Please see the instructions on the software page for setting up a Conda environment based on this file.

Parts of the Notebook use code from a tutorial provided here.

Query data via STAC¶

Import packages.

In [2]:

Copied!





import pystac_client
import stackstac
import pyproj
import xarray
import rioxarray
import pystac_client
import stackstac
import pyproj
import xarray
import rioxarray

Define a point of interest to intersect our query.

In [4]:

Copied!

lon, lat = 11.39, 47.26     # Innsbruck (Austria)
lon, lat = 11.39, 47.26     # Innsbruck (Austria)

Use the pystac_client package and element84’s search endpoint to look for items from the sentinel-s2-l2a-cogs collection on Amazon Web Service (AWS). As search criteria we use a maximum cloud cover (per scene) and a time period.

In [5]:

Copied!

URL = "https://earth-search.aws.element84.com/v0"
catalog = pystac_client.Client.open(URL)
URL = "https://earth-search.aws.element84.com/v0"
catalog = pystac_client.Client.open(URL)

In [11]:

Copied!





%%time                                  # See how fast the query is
items = catalog.search(
    intersects=dict(type="Point", coordinates=[lon, lat]),
    collections=["sentinel-s2-l2a-cogs"],
    query =  ['eo:cloud_cover<50'],
    datetime="2021-04-01/2021-08-31"
).get_all_items()
print(f'{len(items)} items found')      # How many items in the catalog match these criteria?
%%time                                  # See how fast the query is
items = catalog.search(
    intersects=dict(type="Point", coordinates=[lon, lat]),
    collections=["sentinel-s2-l2a-cogs"],
    query =  ['eo:cloud_cover<50'],
    datetime="2021-04-01/2021-08-31"
).get_all_items()
print(f'{len(items)} items found')      # How many items in the catalog match these criteria?

35 items found
CPU times: total: 78.1 ms
Wall time: 1.27 s

Use stackstac to turn those STAC items into a lazy xarray. Without specifying anything, this makes use of dask and reads the data in chunks, with one chunk shaped (1, 1, 1024, 1024) in the time, band, y, x dimensions (Read more about chunk sizes here, here, and here.). Using all the defaults, our data will be in its native coordinate reference system, at the finest resolution of all the assets.

In [12]:

Copied!

%time stack = stackstac.stack(items)
%time stack = stackstac.stack(items)

CPU times: total: 78.1 ms
Wall time: 55.8 ms

Now print info about this xarray.DataArray called 'stack'. Note that everything has run quickly so far, despite the large data volume covered by our query. Up to now, we are only working with metadata, no download or processing of any pixel values.

In [13]:

Copied!

stack
stack

Out[13]:

<xarray.DataArray 'stackstac-f79fa3e66612fef1b91622a0564184e4' (time: 35,
                                                                band: 17,
                                                                y: 10980,
                                                                x: 10980)>
dask.array<fetch_raster_window, shape=(35, 17, 10980, 10980), dtype=float64, chunksize=(1, 1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/26)
  * time                        (time) datetime64[ns] 2021-04-01T10:27:39 ......
    id                          (time) <U24 'S2B_32TPT_20210401_0_L2A' ... 'S...
  * band                        (band) <U8 'overview' 'visual' ... 'WVP' 'SCL'
  * x                           (x) float64 6e+05 6e+05 ... 7.098e+05 7.098e+05
  * y                           (y) float64 5.3e+06 5.3e+06 ... 5.19e+06
    sentinel:latitude_band      <U1 'T'
    ...                          ...
    constellation               <U10 'sentinel-2'
    title                       (band) <U31 'True color image' ... 'Scene Cla...
    common_name                 (band) object None None 'coastal' ... None None
    center_wavelength           (band) object None None 0.4439 ... None None
    full_width_half_max         (band) object None None 0.027 ... None None None
    epsg                        int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

We are only interested in a relatively small area around our main point of interest (defined above for intersection). So we convert the coordinates of this point into the native CRS of the satellite data, add a buffer and then use this rectangle to slice ("clip") the stacked array spatially.

In [14]:

Copied!

x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)     # Re-project the geographic point coordinates to the sat data CRS (i.e. EPSG:32632)
buffer = 5000  # meters
x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)     # Re-project the geographic point coordinates to the sat data CRS (i.e. EPSG:32632)
buffer = 5000  # meters

Now we have a much smaller DataArray (with 1000 by 1000 cells in x/y) called 'stack_aoi'.

In [15]:

Copied!

stack_aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]
stack_aoi
stack_aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]
stack_aoi

Out[15]:

<xarray.DataArray 'stackstac-f79fa3e66612fef1b91622a0564184e4' (time: 35,
                                                                band: 17,
                                                                y: 1000, x: 1000)>
dask.array<getitem, shape=(35, 17, 1000, 1000), dtype=float64, chunksize=(1, 1, 678, 610), chunktype=numpy.ndarray>
Coordinates: (12/26)
  * time                        (time) datetime64[ns] 2021-04-01T10:27:39 ......
    id                          (time) <U24 'S2B_32TPT_20210401_0_L2A' ... 'S...
  * band                        (band) <U8 'overview' 'visual' ... 'WVP' 'SCL'
  * x                           (x) float64 6.758e+05 6.758e+05 ... 6.858e+05
  * y                           (y) float64 5.242e+06 5.242e+06 ... 5.232e+06
    sentinel:latitude_band      <U1 'T'
    ...                          ...
    constellation               <U10 'sentinel-2'
    title                       (band) <U31 'True color image' ... 'Scene Cla...
    common_name                 (band) object None None 'coastal' ... None None
    center_wavelength           (band) object None None 0.4439 ... None None
    full_width_half_max         (band) object None None 0.027 ... None None None
    epsg                        int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

We want only the four 10-m bands and the Scene Classification (SCL) layer.

In [16]:

Copied!

stack_aoi = stack_aoi.sel(band = ['B02', 'B03', 'B04', 'B08', 'SCL'])
stack_aoi
stack_aoi = stack_aoi.sel(band = ['B02', 'B03', 'B04', 'B08', 'SCL'])
stack_aoi

Out[16]:

<xarray.DataArray 'stackstac-f79fa3e66612fef1b91622a0564184e4' (time: 35,
                                                                band: 5,
                                                                y: 1000, x: 1000)>
dask.array<getitem, shape=(35, 5, 1000, 1000), dtype=float64, chunksize=(1, 1, 678, 610), chunktype=numpy.ndarray>
Coordinates: (12/26)
  * time                        (time) datetime64[ns] 2021-04-01T10:27:39 ......
    id                          (time) <U24 'S2B_32TPT_20210401_0_L2A' ... 'S...
  * band                        (band) <U8 'B02' 'B03' 'B04' 'B08' 'SCL'
  * x                           (x) float64 6.758e+05 6.758e+05 ... 6.858e+05
  * y                           (y) float64 5.242e+06 5.242e+06 ... 5.232e+06
    sentinel:latitude_band      <U1 'T'
    ...                          ...
    constellation               <U10 'sentinel-2'
    title                       (band) <U31 'Band 2 (blue)' ... 'Scene Classi...
    common_name                 (band) object 'blue' 'green' 'red' 'nir' None
    center_wavelength           (band) object 0.4966 0.56 0.6645 0.8351 None
    full_width_half_max         (band) object 0.098 0.045 0.038 0.145 None
    epsg                        int32 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5190240.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5300040.00|\n| 0.0...
    resolution:  10.0

Now we trigger the stacking of these selected bands for our area-of-interest by calling compute(), i.e. we load the data into memory. This may take some minutes but then enables us to analyse and plot the data faster. Alternatively (also for larger datasets), we could continue with the lazy array and in subsequent steps only load what is needed.

In [20]:

Copied!

%%time
data = stack_aoi.compute()
%%time
data = stack_aoi.compute()

CPU times: total: 55 s
Wall time: 5min

Try plotting some data (July 2021 only):

In [22]:

Copied!

data.sel(band='B08').sel(time='2021-07').plot.imshow(col='time')
data.sel(band='B08').sel(time='2021-07').plot.imshow(col='time')

Out[22]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4dab57eb0>

No description has been provided for this image

Cloud masking¶

We use the scene classification provided by ESA via the SCL band as a simple (but not the best) means to mask clouds and other problematic pixels.

In [23]:

Copied!





import dask.array as da

scl = data.sel(band=["SCL"])
# Sentinel-2 Scene Classification Map: nodata, saturated/defective, dark, cloud shadow, cloud med. prob., cloud high prob., cirrus
invalid_data = da.isin(scl, [0, 1, 2, 3, 8, 9, 10])
valid_data = data.where(~invalid_data)
# Show some of the cloud masked data
valid_data.sel(band='B08').sel(time='2021-07').plot.imshow(col='time', robust=True)
import dask.array as da

scl = data.sel(band=["SCL"])
# Sentinel-2 Scene Classification Map: nodata, saturated/defective, dark, cloud shadow, cloud med. prob., cloud high prob., cirrus
invalid_data = da.isin(scl, [0, 1, 2, 3, 8, 9, 10])
valid_data = data.where(~invalid_data)
# Show some of the cloud masked data
valid_data.sel(band='B08').sel(time='2021-07').plot.imshow(col='time', robust=True)

Out[23]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4d9e083d0>

Let's show cloud-masked RGB true color composites.

In [24]:

Copied!

valid_data.sel(time='2021-07').sel(band=['B04', 'B03', 'B02']).plot.imshow(col='time', robust=True)
valid_data.sel(time='2021-07').sel(band=['B04', 'B03', 'B02']).plot.imshow(col='time', robust=True)

Out[24]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4d9e133a0>

Spectral index calculation (NDVI)¶

We calculate the Normalized Difference Vegetation Index (NDVI). If you need more indices consider the Spyndex package.

In [25]:

Copied!

# Compute the NDVI from the NIR (B8) and red (B4) bands
nir, red = valid_data.sel(band='B08'), valid_data.sel(band='B04')
NDVI = (nir - red) / (nir + red)
# Compute the NDVI from the NIR (B8) and red (B4) bands
nir, red = valid_data.sel(band='B08'), valid_data.sel(band='B04')
NDVI = (nir - red) / (nir + red)

NDVI is a new DataArray. If you want to select the NDVI raster closest to a specific date and export it as a GeoTIFF, just run this line:

In [26]:

Copied!

NDVI.sel(time='2021-07-15', method='nearest').rio.to_raster("C:\work\etrainee\stac_exports/S2-NDVI_midjuly.tif")
NDVI.sel(time='2021-07-15', method='nearest').rio.to_raster("C:\work\etrainee\stac_exports/S2-NDVI_midjuly.tif")

Let's plot one month of NDVI rasters.

In [27]:

Copied!

NDVI.sel(time='2021-07').plot.imshow(col='time', cmap="PRGn", robust=True)
NDVI.sel(time='2021-07').plot.imshow(col='time', cmap="PRGn", robust=True)

Out[27]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4da128bb0>

Time series operations¶

Interpolate missing values¶

Let's fill missing values by (linear) interpolation between observations.

In [28]:

Copied!

NDVI.sel(time='2021-07').interpolate_na(dim='time').plot(col='time', cmap="PRGn", robust=True)
NDVI.sel(time='2021-07').interpolate_na(dim='time').plot(col='time', cmap="PRGn", robust=True)

Out[28]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4da070f40>

Resampling¶

We resample/aggregate to one month temporal resolution by computing the mean for each month. Other aggregations (e.g. maximum) can be achieved similarly. Note how we can absolutely define the range of values to be displayed (important for visual comparisons).

In [35]:

Copied!

NDVI.resample(time='1M').mean().plot(col='time', cmap="PRGn", vmin=-0.2, vmax=0.8)
NDVI.resample(time='1M').mean().plot(col='time', cmap="PRGn", vmin=-0.2, vmax=0.8)

Out[35]:

<xarray.plot.facetgrid.FacetGrid at 0x1b4923307f0>

Trend analysis¶

The xarrayutils package contains a convenient linear_trend() function. We use this to fit a linear least-squares trendline to our NDVI time series.

In [38]:

Copied!

from xarrayutils import linear_trend
NDVI_regressed = linear_trend(NDVI, 'time')
from xarrayutils import linear_trend
NDVI_regressed = linear_trend(NDVI, 'time')

Plot maps (with default colormaps) for the resulting slope of the estimated linear trend along with r-value (Pearson correlation coefficient), p-value (for a hypothesis test whose null hypothesis is that the slope is zero), and standard error (of the estimated slope). We use matplotlib syntax and functions to configure the plot.

In [50]:

Copied!





from matplotlib import pyplot as plt
fig, axes = plt.subplots(2,2, figsize=(15,10), sharex=True, sharey=True)
NDVI_regressed.slope.plot(robust=True, ax=axes[0,0])
NDVI_regressed.r_value.plot(robust=True, ax=axes[0,1])
NDVI_regressed.p_value.plot(robust=True, ax=axes[1,0])
NDVI_regressed.std_err.plot(robust=True, ax=axes[1,1])
for ax, param in zip(axes.flat, ['slope', 'r-value', 'p-value', 'std error']):
    ax.set_title(f'NDVI trend 04-08/2021: {param}')
fig.tight_layout()
from matplotlib import pyplot as plt
fig, axes = plt.subplots(2,2, figsize=(15,10), sharex=True, sharey=True)
NDVI_regressed.slope.plot(robust=True, ax=axes[0,0])
NDVI_regressed.r_value.plot(robust=True, ax=axes[0,1])
NDVI_regressed.p_value.plot(robust=True, ax=axes[1,0])
NDVI_regressed.std_err.plot(robust=True, ax=axes[1,1])
for ax, param in zip(axes.flat, ['slope', 'r-value', 'p-value', 'std error']):
    ax.set_title(f'NDVI trend 04-08/2021: {param}')
fig.tight_layout()

Create an interactive map to inspect the slope of the trend in more detail.

In [ ]:

Copied!

import hvplot.xarray
import holoviews as hv
NDVI_regressed.slope.hvplot(x='x', y='y', aspect=1, clim=(-0.02,0.02), cmap="PRGn", title='NDVI slope April - Aug')
import hvplot.xarray
import holoviews as hv
NDVI_regressed.slope.hvplot(x='x', y='y', aspect=1, clim=(-0.02,0.02), cmap="PRGn", title='NDVI slope April - Aug')

Create an animation¶

Let's create an animation of our NDVI time series using the xmovie package. We resample to weekly means of NDVI and interpolate missing values.

In [58]:

Copied!





%%time
NDVI_weekly = NDVI.resample(time='1W').mean()
NDVI_weekly = NDVI_weekly.chunk(dict(time=-1))                        # rechunk into a single dask array chunk along the time dimension
NDVI_weekly = NDVI_weekly.interpolate_na(dim='time')
%%time
NDVI_weekly = NDVI.resample(time='1W').mean()
NDVI_weekly = NDVI_weekly.chunk(dict(time=-1))                        # rechunk into a single dask array chunk along the time dimension
NDVI_weekly = NDVI_weekly.interpolate_na(dim='time')

CPU times: total: 141 ms
Wall time: 149 ms

Now be patient while the movie is created.

In [59]:

Copied!

%%time
from xmovie import Movie
mov = Movie(NDVI_weekly.chunk({'time':1}), dpi=100, vmin=-0.8, vmax=0.8, cmap="PRGn")
%%time
from xmovie import Movie
mov = Movie(NDVI_weekly.chunk({'time':1}), dpi=100, vmin=-0.8, vmax=0.8, cmap="PRGn")

CPU times: total: 1min 15s
Wall time: 1min 2s

Use parallelization by Dask for saving the frames to an MP4 video clip (or a GIF). For this to work we have chunked the input DataArray with a single step along the dimension that represent our frames. Again be patient ... (Code cell doesn't stop running even after video clip is saved!?)

In [ ]:

Copied!

%%time
mov.save(r'C:\work\etrainee\stac_exports\NDVI_animation.mp4', overwrite_existing=True, parallel=True, framerate=3)
%%time
mov.save(r'C:\work\etrainee\stac_exports\NDVI_animation.mp4', overwrite_existing=True, parallel=True, framerate=3)